#include "gctl/core.h"
#include "gctl/io.h"

#include "../utilities/efem_sf.h"

using namespace gctl;

void get_corner(std::string poly_file, array<point2dc> &cor)
{
	int tmp_id;
	std::ifstream infile(poly_file);

	infile >> tmp_id >> tmp_id >> tmp_id >> tmp_id;

	cor.resize(3, point2dc(0.0, 0.0));
	for (int i = 0; i < 3; ++i)
	{
		infile >> tmp_id >> cor[i].x >> cor[i].y;
	}

	infile.close();
}

int main(int argc, char const *argv[]) try
{
	/*
	int edge_id = atoi(argv[1]);
	std::string tri_file = "data/tri.1";
	std::string poly_file = "data/tri.poly";

	array<vertex2dc> nodes;
	array<triangle2d> eles;
	gctl::read_Triangle_node(tri_file, nodes);
	gctl::read_Triangle_element(tri_file, eles, nodes);

	array<point2dc> cor;
	get_corner(poly_file, cor);

	array<point3dc> nodeval(nodes.size(), point3dc(0.0,0.0,0.0));
	array<double> nodex(nodes.size());
	array<double> nodey(nodes.size());
	array<double> edge_tang(nodes.size(), 0.0);

	gmshio fio;
	fio.init_file(tri_file, Output);
	fio.set_packed(NotPacked, Output);
	fio.save_mesh(eles, nodes);

	double edge_x, edge_y;
	int tmp_id;
	if (edge_id <= 5) tmp_id = edge_id/2;
	else if (edge_id == 6) tmp_id = 1;
	else if (edge_id == 7) tmp_id = 2;

	edge_x = cor[(tmp_id+1)%3].x - cor[tmp_id].x;
	edge_y = cor[(tmp_id+1)%3].y - cor[tmp_id].y;

	double evax, evay;
	//efem_sf::linear elsf;
	efem_sf::quadratic eqsf;
	for (int i = 0; i < nodes.size(); ++i)
	{
		//elsf.triangle(nodes[i].x, nodes[i].y, cor[0].x, cor[1].x, cor[2].x, cor[0].y, cor[1].y, cor[2].y, 
		//	edge_id, efem_sf::Value, efem_sf::Normal, evax, evay);

		eqsf.triangle(nodes[i].x, nodes[i].y, cor[0].x, cor[1].x, cor[2].x, cor[0].y, cor[1].y, cor[2].y, 
			edge_id, efem_sf::Value, efem_sf::Normal, evax, evay);

		nodeval[i].x = evax;
		nodeval[i].y = evay;

		nodex[i] = evax;
		nodey[i] = evay;

		edge_tang[i] = (evax*edge_x + evay*edge_y)/sqrt(edge_x*edge_x + edge_y*edge_y);
	}
	fio.save_data("Node Vector", nodeval, NodeData);
	fio.save_data("Node X", nodex, NodeData);
	fio.save_data("Node Y", nodey, NodeData);
	fio.save_data("Edge tang", edge_tang, NodeData);
	*/

	
	double Ex = 2.0, Ey = 1.2;

	std::string tri_file = "data/tri.1";
	std::string poly_file = "data/tri.poly";

	array<vertex2dc> nodes;
	array<triangle2d> eles;
	gctl::read_Triangle_node(tri_file, nodes);
	gctl::read_Triangle_element(tri_file, eles, nodes);

	array<point2dc> cor;
	get_corner(poly_file, cor);

	double edge_E[8];
	double edge_x, edge_y;
	for (int i = 0; i < 3; ++i)
	{
		edge_x = cor[(i+1)%3].x - cor[i].x;
		edge_y = cor[(i+1)%3].y - cor[i].y;
		edge_E[2*i] = edge_E[2*i+1] = (edge_x*Ex + edge_y*Ey)/sqrt(edge_x*edge_x + edge_y*edge_y);
	}
	edge_E[6] = edge_E[2];
	edge_E[7] = edge_E[4];

	array<point3dc> nodeval(nodes.size(), point3dc(0.0,0.0,0.0));
	array<double> nodex(nodes.size(), 0.0);
	array<double> nodey(nodes.size(), 0.0);

	gmshio fio;
	fio.init_file(tri_file, Output);
	fio.set_packed(NotPacked, Output);
	fio.save_mesh(eles, nodes);

	double evax, evay;

	//efem_sf::linear elsf;
	//for (int i = 0; i < nodes.size(); ++i)
	//{
	//	for (int j = 0; j < 3; ++j)
	//	{
	//		elsf.triangle(nodes[i].x, nodes[i].y, cor[0].x, cor[1].x, cor[2].x, cor[0].y, cor[1].y, cor[2].y, 
	//			j, efem_sf::Value, efem_sf::Normal, evax, evay);
//
	//		nodeval[i].x += evax*edge_E[j];
	//		nodeval[i].y += evay*edge_E[j];
//
	//		nodex[i] += evax*edge_E[j];
	//		nodey[i] += evay*edge_E[j];
	//	}
	//}

	efem_sf::quadratic eqsf;
	for (int i = 0; i < nodes.size(); ++i)
	{
		for (int j = 0; j < 8; ++j)
		{
			eqsf.triangle(nodes[i].x, nodes[i].y, cor[0].x, cor[1].x, cor[2].x, cor[0].y, cor[1].y, cor[2].y, 
				j, efem_sf::Value, efem_sf::Normal, evax, evay);

			nodeval[i].x += evax*edge_E[j];
			nodeval[i].y += evay*edge_E[j];

			nodex[i] += evax*edge_E[j];
			nodey[i] += evay*edge_E[j];
		}
	}

	fio.save_data("Node Vector", nodeval, NodeData);
	fio.save_data("Node X", nodex, NodeData);
	fio.save_data("Node Y", nodey, NodeData);
	
	return 0;
}
catch (std::exception &e)
{
	std::cerr << e.what() << std::endl;
}